AB-INITIO GUTZWILLER METHOD: FIRST APPLICATION 
TO PLUTONIUM 



J.-P. JULIEN 

Laboratoire d'Etudes des Proprietes Electroniques des Solides, 
CNRS, P.P. 166, 38042 Grenoble Cedex 9, France 

AND 

J. BOUCHET 
CEA-DAM, DPTA, 
Bruyeres-le-Chatel, France 



Abstract. Using a density matrix approach to Gutzwiller method, we 
present a formalism to treat ab-initio multiband Tight-Binding Hamilto- 
nians including local Coulomb interaction in a solid, like, for e.g., the de- 
generate Hubbard model. We first derive the main results of our method: 
starting from the density matrix of the non-interacting state, we build a 
multi-configurational variational wave function. The probabilities of atomic 
configurations are the variational parameters of the method. The kinetic 
energy contributions are renormalized whereas the interaction contribu- 
tions are exactly calculated. A renormalization of effective on-site levels, 
in contrast to the usual one-band Gutzwiller approach, is derived. After 
minimization with respect to the variational parameters, the approximate 
ground state is obtained, providing the equilibrium properties of a material. 
Academic models will illustrate the key points of our approach. Finally, as 
this method is not restricted to parametrized Tight-Binding Hamiltonians, 
it can be performed from first principles level by the use of the so-called 
"Linearized Muffin Tin Orbitals" technique. To avoid double counting of the 
repulsion, one subtracts the average interaction, already taken into account 
in this density functional theory within local density approximation (DFT- 
LDA) based band structures method and one adds an interaction part " a la 
Hubbard" . Our method can be seen as an improvement of the more popular 
LDA+C7 method as the density-density correlations are treated beyond a 
standard mean field approach. First application to Plutonium will be pre- 
sented with peculiar attention to the equilibrium volume, and investigations 
for other densities will be discussed. 
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1. Introduction 

Except for small molecules, it is impossible to solve many electrons sys- 
tems without imposing severe approximations. If the configuration interac- 
tion approaches (CI) or Coupled Clusters techniques [1] are applicable for 
molecules, their generalization for solids is difficult. For materials with a 
kinetic energy greater than the Coulomb interaction, calculations based on 
the density functional theory (DFT), associated with the local density ap- 
proximation (LDA) [2,3] give satisfying qualitative and quantitative results 
to describe ground state properties. These solids have weakly correlated 
electrons presenting extended states, like sp materials or covalent solids. 
The application of this approximation to systems where the wave func- 
tions are more localized (d or /-states) as transition metals oxides, heavy 
fermions, rare earths or actinides is more questionable and can even lead to 
unphysical results : for example, insulating FeO and CoO are predicted to 
be metalic by the DFT-LDA. On another hand, theoretical "many body" 
approaches like diagrammatic developments [4], slave bosons [5], decoupling 
of the equations of motion of Green functions by projection techniques [6], 
and more recently dynamical mean field theory (DMFT) in infinite dimen- 
sion [7], treat in a much better way correlation effects than the DFT-LDA 
does. However the price to be paid is an oversimplification of the system, 
generally reducing the number of involved orbitals and using parameterized 
Hamiltonians (like Hubbard model) where the ab-initio aspect of the DFT- 
LDA is lost. Finally, these methods, contrary to the DFT-LDA, are scarcely 
variational. Recently several attemps have been proposed to couple these 
two points of view as in the LDA+U [8], or LDA+DMFT [9,10] approaches. 
In the same spirit, the approach we describe below, tries to keep advantages 
on both sides: it is a variational method which is multi-configurational, 
contrary to the DFT-LDA, but without loosing the "adjustable parameters 
free" advantages of the ab-initio side. The next section is devoted to the 
derivation of our formalism. It is then applied to known academic cases 
to prove the reliability of our approach. The insertion of this approach at 
the ab-initio level is presented in section 2.5. The nature of the electronic 
structure of Plutonium being still under discussion, the application of our 
method, in the last section, is in accordance with previous works and also 
gives some new insights for this material. 

2. Method 

2.1. GUTZWILLER APPROACH FOR THE ONE-BAND HUBBARD MODEL 

Among numerous theoretical approaches, the Gutzwiller method [11, 12] 
provides a transparent physical interpretation in term of atomic configu- 
rations of a given site. Originally it was applied to the one-band Hubbard 
model Hamiltonian [13]: 
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(1) 
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(2) 



and 



U n i\ n ii 



(3) 



which contains a kinetic part H^m with a hopping integral Uj from site 
j to site i, and an interaction part with a local Coulomb repulsion U for 
electrons on the same site. c] CT (respectively Cj CT ) is the creation (respectively 
destruction) operator of an electron at site i with up or down spin a. = 
cj^Cicr measures the number (0 or 1) of electron at site i with spin a. This 
Hamiltonian contains the key ingredients for correlated up and down spin 
electrons on a lattice: the competition between derealization of electrons by 
hoppings and their localization by the interaction. It is one of the most used 
models to study electronic correlations in solids (for a review see Ref. [14]). 

In the absence of the interaction U, the ground state is that of uncor- 
rected electrons |^o) an d has the form of a Slater determinant. As U is 
turned on, the weight of doubly occupied sites must be reduced because 
they cost an additional energy U per site. Accordingly, the trial Gutzwiller 
wave function (GWF) \^g) is built from the Hartree-like uncorraleted wave 
function (HWF) |* >, 



The role of g is to reduce the weight of configurations (i.e. a way of 
spreading N electrons over the lattice) with doubly occupied sites, where 
D = J2i n i\ n i\ measures the number of double occupations and g (< 1) is a 
variational parameter. In fact, this method corrects the mean field (Hartree) 
approach for which up and down spin electrons are independent, and, some 
how, overestimates configurations with double occupied sites. Using the 
Rayleigh-Ritz principle, this parameter is determined by minimization of 
the energy in the Gutzwiller state |^g), giving an upper bound to the true 
unknown ground state energy of H. Note that to enable this calculation to 
be tractable, it is necessary to use the Gutzwiller's approximation which 
assumes that all configurations in the HWF have the same weight. Details 
of the derivation can be found in the article of Vollhardt [15]. 

Nozieres [16] proposed an alternative way, showing that the Gutzwiller 
approach is equivalent to renormalize the density matrix in the GWF which 
can be reformulated as: 



\v G ) = g D \*o) 



(4) 
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PG 



T ] Po T (5) 



The density matrices pq = \^g){' i Pg\ an d po = |V'o)(V ; o| are projectors 
on the GWF and HWF respectively. T is an operator, diagonal in the 
configuration basis, T = FJi T{ and Tj is a diagonal operator acting on site 

i 

T i \L i ,L') = J^f\L h L') (6) 
V Po{Li) 

Li is an atomic configuration of the site i, with probability p(Li) in the 
GWF and po(Li) in the HWF respectively, whereas L' is a configuration 
of the remaining sites of the lattice. Note that this prescription does not 
change the phase of the wave function as the eigenvalues of the operators 
Tj are real. The correlations are local, and the configuration probabilities 
for different sites are independent. 
The expectation value of Hamiltonian (1) is given by 

(H) G = Tr( PG H) (7) 

The mean value of one-site operators (interaction U) is exactly calculated 
with the double occupancy probability di = (n^rii^G- di is the new vari- 
ational parameter replacing g. From expressions (5) and (6), the two-sites 
operators contributions of the kinetic energy can be written as 



(cU 3a ) G = Tr( PG ^) = (4^)0 E J J (8) 

where L' a and L a are the only two configurations of spin a at sites i and 
j that give non-zero matrix element to the operator in the brackets as 
illustrated on Fig.l. The summation is performed over the configurations of 
opposite spin L_ CT . Their corresponding probabilities are pictured on Table 
1 for an homogeneous state (for any site i, {rii a ) = n and {n^riii) = d). 
The probabilities po in the HWF depend only on the number of electrons, 
whereas the p in the GWF also depend on di. 

After some elementary algebra, one can show that the Gutzwiller mean 
value can be factorised: 

( C l C ^>G = V / ^(4a C i^)0 V / ^ : ( 9 ) 

where these renormalization factors qi a are local and can be expressed as: 
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Figure 1. Initial and final configurations contributing to the mean value {c ia Cj a )G- 



Li po(Li) 



p(Li) 



(1-n) 2 l-2n + d 



| n(l — n) n — d 



I n(l-n) 



-d 



TABLE 1. Different possi- 
ble configurations of one site 
and the corresponding proba- 
bilities in the HWF (po) and 
GWF (p). 



where Tli(J IS Hi shorthand for (nj CT ), the average number of electron on the 
considered "orbital-spin" in the HWF, which could be site and/or spin 
independent if the state is homogeneous and/or paramagnetic (it is the case 
we consider here for pedagogy, dropping indices ia). The kinetic energy e^ in 
of the non-interacting electrons state is renormalized by a factor q which 
is smaller than one in the correlated state, and equal to one in the HWF. 
Then, we minimize the variational energy 



6 



J.-P. JULIEN AND J. BOUCHET 



E(d) = (H) G = qe° km + Ud (11) 

with respect to d. In the case of half filling (n = 1/2), if the repulsion U 
exceeds a critical value U c = 8e^ n , q is equal to zero, leading to an infinite 
quasiparticle mass with a Mott-Hubbard Metal-Insulator transition which 
is, in this context, often referred to as "the Brinkmann-Rice transition" [17], 
as these authors first applied the Gutzwiller approximation to the Metal- 
Insulator transition. Application of this "one orbital per site" formalism for 
inhomogeneous states is possible because all involved quantities are local. 
An example can be found in [18] for model CuOi planes, in connection 
with the electronic structure of High Tq superconductors. 

2.2. INEQUIVALENT SITES: RENORMALIZATION OF LEVELS 

When sites are inequivalent, or if orbitals belong to different symmetries as 
in a multiorbital spdf basis case of further sections, it is necessary to add 
to the Hamiltonian an on-site energy term 

H n — site = ^ ] ^ia^icr (12) 
icr 

Hence this enlarged Hubbard Hamiltonian can be written as 

H = J2 *ij4r c Jff + S + U J2 n H n U ( 13 ) 

i^j : cj i<7 i 

In that case, the starting HWF, directly obtained from the non-interacting 
part of the Hamiltonian, is not automatically the best choice, giving the 
optimal GWF, i.e. having the lowest energy. For example, if we look for the 
ground state of Hamiltonian (13) in the Hartree-Fock (HF) self-consistent 
field formalism, it is necessary to vary the orbital occupations. Practically, it 
can be achieved by replacing this Hamiltonian, by an effective Hamiltonian 
H e f f of independent particles with renormalized on-site energies : 

R eff = Yl UAcja + ^n i(T (+C) (14) 
i^j,a ia 

The HWF we are looking for, is an approximate ground state of the 
true many-body Hamiltonian (13) and is the exact ground state of effective 
Hamiltonian (14). The additive constant C accounts for double counting 
energy reference, so that the ground state energies are the same for both 
Hamiltonians: 

(H ef f) = (H) (15) 

The effective Hamiltonian depends on pararemeters ej CT . The optimal choice 
can be obtained by minimizing the ground state energy of H e ff with respect 
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to these parameters. With the help of Hellmann-Feyman theorem, one can 
easily see that the derivative of the kinetic energy is 

djHkin) _ \- d(n ja ) 

Be ~ ^ ja Be ( ' 

&i,o 061(7 

On another hand, differentiation of equality (15) associated with expression 
(16) and the mean field approximation (n^n^) rs ( n it)( n i|) enables to 
retrieve the well-known formula for the on-site energies 



e ia = e° ia + U{m- a ) (17) 

and the constant C is simply — U Y,i{ n il) i n il) ■ 

In the Gutwiller approach, the same argument about the variation of 
orbital occupation, i.e. flexibility on the HWF |\l/o)> is true. It is necessary 
to find a way to vary this Slater determinant from which the GWF |*g) is 
generated, so that the Gutzwiller ground-state energy is minimum. Clearly 
one has to find an equivalent of formula (17) in the Gutzwiller context, 
which has never been established, to our knowledge. The average value of 
Hamiltonian (13) on a GWF is given by: 



(* G \H\* G ) = £ 4^(4^)0^ + ^ £ <*» + £ 4Kr>o (18) 

ijcr i iacr 

Following the same footing of previous HF self-consistent field approach, 
one has to find an effective Hamiltonian H e ff of independent particles hav- 
ing |^o) as exact ground state. This state |^o) generates the GWF 
which is an approximate ground state of the true interacting Hamiltonian 
(13). In analogy with (15), the condition 

<tf |#e//|*0> = (Vg\H\V G ) (19) 

leads to the expression for the searched H e ff. 

H eff = £ UAcjc + £ ^iaUia + C (20) 
iy^j,cr icr 

with effective but fixed renormalized hoppings Uj = ^fq^^ij ^Qja and having 
effective on-site energies ej CT which have still to be determined. Hellmann- 
Feynman theorem applied to H e ff provides again an expression similar to 
(16), but with effective hoppings. Taking into account the dependence of 
the qia's through nj CT in differentiating (18) and (19) with respect to the pa- 
rameters ei a , after some calculations, one obtains the equivalent expression 
of (17) in the Gutzwiller context: 
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dn ia 



(21) 



Here ej CT is the partial kinetic energy of orbital-spin ia, it is given by 



with N ia the icr-projected density of states (DOS) for Hamiltonian H e ff. 
The remaining constant C' that ensures (19) explicitly reads 



To solve the full problem of finding an approximate ground state to 
Hamiltonian (13), one is faced to a self-consistent loop which can be pro- 
ceeded in two steps. First one can get the occupations (nj CT )o from a HWF, 
and a set of 'bare' e® a levels. Then one obtains a set of configuration pa- 
rameters, the probabilities of double occupation, di by minimizing (18) with 
respect to these probabilities. Afterwards the on-site levels are renormal- 
ized according to (21) and the next loop starts again for the new effective 
Hamiltonian H e ff till convergence is achieved. 

As illustration of the importance of the renormalization of levels, we 
studied the case of an alternate Hubbard chain (often called Ionic Hubbard 
model). This model contains two kinds of alternating atoms A and B on 
an infinite chain, having on-site energies and e° B respectively, coupled 
by a hopping integral t. The same local Coulomb repulsion U acts on each 
site. For a given total number of electrons, one can fix a repartition of 
electrons among sites A and B, and compute the energy of the ground 
state: first within the HF mean field approximation, and secondly, within 
the Gutzwiller approach. Browsing the electronic occupation of A-site, by 
adjunction of Lagrange multiplier to fix it to a given value, one looks for the 
lowest energy state. The corresponding ground state energies as function of 
the A-filling are presented on Fig. 2. First of all, the lowest Hartree state 
could be more efficiently directly found, after some self-consistent loops, 
via the on-site renormalization of levels of equation (17) in the HF context, 
as explained in precedent paragraph. By inspection of the curve, it is also 
obvious that this lowest Hartree does not generate the lowest Gutzwiller 
state. It is necessary to browse among different A-fillings to find the best 
Gutzwiller ground state. If this browsing procedure is still tractable for 
simple models, as we did in Ref. [18], its generalization to multiorbitals 
cases would be practically impossible. It is the main advantage of formula 
(21) to avoid this cumbersome search for optimized levels and to provide 
a systematic way of finding them, similar to (17), leading to the best (i.e. 
lowest) Gutzwiller ground state. 




(22) 




(23) 
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Figure 2. Total energy of the alternate chain versus A-subband filling. Upper curve: 
Hartree-Fock result, Lower curve: Gutzwiller result. The 2 minima are clearly different 



2.3. GENERALIZATION TO THE DEGENERATE HUBBARD 
HAMILTONIAN 

Now we generalize this density matrix formalism [19, 20] for the degenerate 
Hubbard Hamiltonian which, with usual notations, reads: 



H = tiajpCiaaCjfa + H int 



(24) 



with the model interaction 



hit 



E 

i,acrj^f3cr' 



(25) 



where a, f3 and a, a' are orbitals and spins index respectively, necessary 
to account for orbital degeneracy. The case ao = (3a' is excluded from the 
interaction because of Pauli principle. We neglect any spin flip term in the 
interaction for simplicity. They could be in principle taken into account in 
our approach, as it is done in a different work by Bunemann et al [21]. 
However this procedure would involve a diagonalization of atomic part of 
Hamiltonian, that complicates the presentation of our approach without 
bringing any new physical ingredients. 
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As in the one-band case, we define a Gutzwiller renormalized density 
matrix with the operator T given by eqs (5) and (6). The main difference 
being the greater number of atomic configurations, equal to 2 2N , with N the 
orbital degeneracy For a given site we have now probabilities for double, 
triple, etc... multiple occupancy, which are the new variational parameters 
generalizing the role of d. Of course, the number of independent variational 
probabilities is smaller than the number of configurations, as different con- 
figurations could have the same probabilities for symmetries reasons. For 
example, in a paramagnetic case, a configuration and its spin reverse are 
equivalent leading to the same probability. Moreover, the probabilities are 
not independent of each other as the sum over all probabilities have to be 
equal to 1, and we have also to conserve the average electronic occupa- 
tion of given orbital-spin (nj aCT ). These constraints could be either directly 
included in the expressions of empty and single occupied configurations 
probabilities, or treated by adjunction of Lagrange multipliers as in the 
slaves bosons approach [5] . This last formulation has the advantage of giv- 
ing more symmetric expressions. Using the expression (6) of Tj operators, 
we can directly obtained the factorized form of the kinetic energy terms: 

( c Lt<Wg = VQi^(4aa C jl3a)0^qj^ (26) 

where the g-factors reduce the kinetic energy and are expressed as functions 
of the variational parameters and the number of electrons according to: 

= === == \Jp{ia<? ■ unocc, L'^)p{iaa : occ, L[) (27) 

i 

Here p(iaa : occ, (respectively p(iaa : unocc, LQ) represents the proba- 
bility of the atomic configuration of site i, where the orbital a with spin a is 
occupied (resp. unoccupied) and where L[ is a configuration of the remain- 
ing orbitals of this site. This result is similar to the expression obtained 
by Biinemann et al. [22], but it is obtained more directly by the density 
matrix renormalization (5). To obtain the expression of the qi aa factors, an 
additional approximation to the density matrix of the uncorrelated state 
was necessary. This approximation can be viewed as the multiband gener- 
alization of the Gutzwiller approximation, exact in infinite dimension [23] 

(LL"\ Po \L'L") « Po (L")J2(LL"\p \L'L") (28) 

L" 

Where we have replaced an off-diagonal element of the density by its aver- 
age value over the configurations L" . L and L' are configurations of one or 
two sites, involved in the calculation of interaction or kinetic term and L" 
is the configuration of remaining sites. This approximation allows to per- 
form calculations, and however preserves sum rules of the density matrix. 
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Similarly, the interaction between an electron at site i on orbital a with 
spin a and an electron on orbital (5 with spin a' involves a term 



where L\ is a configuration of the remaining spin-orbitals of this site, other 
than aa and (5a 1 . 

As illustration, we studied the academic case of paramagnetic state for 
doubly degenerate bands like, for instance, e s -symmetry ci-orbitals in cubic 
or octahedral environment. Hybridization among these degenerate orbitals 
is supposed to produce a kinetic energy e\ in in the uncorrelated state. We 
take a model interaction where the general expression (25) reduces to: 



with two independent parameters U and U' as the relation U — U' = 2 J 
stands [24]. The interaction between electrons of same spin is reduced by 
the exchange integral J, which is essential to reproduce first Hund's law of 
maximum spin. The application of the above prescription directly leads to 
the variational energy: 



E G = 2qe° kin + 2Ud + 2U'd 1 + 2{U' - J)d 2 + 2{U + 2U' - J) (24 + /) (31) 



/ and t are the quadruple and triple occupancy respectively, whereas there 
are three possibilities of double occupancies : do (same orbital, different 
spin), d\ (different orbital, different spin)and d 2 (different orbital, same 
spin). Some of the corresponding configurations with multiple occupancy 
are pictured on Table 2, followed by their probability and their interac- 
tion energy. This expression is identical to the result obtained by different 
authors using Gutzwiller-type wave function [25,26], or multiband slave- 
boson approach [27] which is a multiorbital generalization of Ref. [5]. It is 
to be stressed the very physical "transparent" approach with the density 
matrix formalism, leading to simple expressions. Also, there is no approxi- 
mation about less favorable configurations, discarded from the beginning as 
in Ref. [28] . For a given electronic filling, we use a Newton- Raphson proce- 
dure to minimize Eq with respect to do, d\, d 2 , t and /. We again choose a 
half filled case, and we scale all contributions in term of the kinetic energy. 

On Fig. 3 we plot the probabilities of different configurations versus the 
direct Coulomb interaction U. It can be seen that the system undergoes a 
metal- insulator transition for a sufficiently high value of U, close to 9. It is 
easy to perform the same kind of calculation in the case of triply degenerate 




(29) 



Hint — U 2~2iaa n iatr n ia—cr ~\~ U ^2ia^/3a n iaa n ij3—o 

+(U> - J) 



(30) 
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do U 



Tl 






T 


1 




T 


T 




Tl 


1 




Tl 


Tl 



TABLE 2. 



di U' 



d 2 U' - J 



t U + 2U' — J 



f 2U + 2(2(7' - J) 




Figure 3. For the half filling case all the probabilities are equal for (7 = 0. For U c ~ 9 
we have a transition from the metallic to the insulator state which is of first order at half 
filling and second order for other concentration as seen for n — 0.48. 



orbitals relevant for d-orbitals with the symmetry t% g , (for instance Ti t2 g 
in LaTiOs), /-orbitals with the symmetry T± or T2 in cubic or octahedral 
environment (rare-earth element or actinides) or p-orbitals like in fullerene 

C60- 
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2.4. DIFFERENT ORBITAL SYMMETRIES: FULLY HYBRIDIZED 
HAMILTONIAN 

As our aim is to describe realistic materials from the ab-initio level, it is 
necessary to have a full set of spd and possibly / (for actinides or rare 
earths) basis. The system can be described by the following Hamiltonian 
which is the sum of a kinetic term, local Coulomb repulsions and an on-site 
contribution accounting of difference of site (and/or symmetries) 



The kinetic term is responsible of the hybridization of orbitals with 
different ^-symmetries on neignboring sites among each others. Moreover, we 
assume that the interaction part of the Hamiltonian Hi n t only concerns 
one subset of correlated orbitals (say /). All atomic configurations T of 
this subset on each site i will be considered. Note that H can be seen as 
a multiband hybrid between the Hubbard Hamiltonian and the periodic 
Anderson Hamiltonian. It contains hybridization of localized interacting / 
orbitals among each others (Hubbard) but also with extended spd states 
(Anderson). Using the results of previous sections, the variational energy 
in the Gutzwiller state reads 



In this expression the q- factors of site i are functions, through (27), of the 
probabilities Pi(T) of the atomic configurations T of /-orbitals at the same 
site. They are equal to 1 if orbital a or (5 does not belong to this subset, 
i.e. for extended states. U-p is a proper combination of Coulomb direct and 
exchange contributions U aa p a i , accounting for the interaction energy which 
arises as prefactors of expressions (29), and which can be seen for e.g. on 
simplified case of (31). As in our previous simpler models, the probabilities 
Pi(T) are the variational parameters and one has to minimize Eq with 
respect to each of them (and at each inequivalent site) according to 



H = ^li^ja(3iT tiot,j/3 c i aa c jl3(T 



(32) 



iota ^iaa( n icu^)o 



(33) 





l-(\(7 

d Pi (T) 



(34) 
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To avoid the cumbersome calculations of all (c] acr Cj i a cr )o cross-terms which 
are present in a fully hybridized case, we propose a recursive procedure for 
the minimization of Eq, i.e. the search for the optimal set of probabilities. In 
spirit close to the impurity model of the DMFT approach, we first consider 
that there is only one correlated site, say i = 0. This site is supposed 
to be embedded in a reference fixed medium where all q's other than the 
considered site '0', are equal to 1 at beginning or to their previous values in 
the self-consistent process that has to be performed afterwards. Then the 
set of equations (34) for all configurations Tq of site can be rewritten as 



dE G _ ^ dln{y/%. 



The partial kinetic energies ei = o aa of orbital aa at site i = 0, gener- 
alizing (22), are obtained from partial projected DOS, available from any 
electronic structure code, and computed with site i = embedded in the 
reference medium. To ensure numerical stability when solving system of 
eqs. (35) and according to the spirit of Landau theory of Fermi liquids, 
the interactions are progressively switched on from zero to their final val- 
ues starting from the probabilities of uncorrelated case. After solution, the 
probabilities Pq{T) are used to compute the local g-factors of site 0. If all 
sites are equivalent, one would get the same results on other sites. Accord- 
ingly, the g-factors of other sites are all set equal to the 'O'-th ones (one 
would have to repeat this impurity-like calculation if there are inequivalent 
sites, i.e. crystal structures with more than one atom per cell or disordered 
systems). Changing the g-factors affects the partial kinetic energies ej=oa<j 
and also the occupation of orbitals, as the reference medium now has new 
effective hoppings. Process must be iterated till convergence. The advan- 
tage of this way of solving iteratively eqs. (34) is that the only required 
ingredients to get the solutions are partial (local) kinetic energies and oc- 
cupations of orbitals at site directly obtained from partial DOS's. The 
price to be paid is a greater number of electronic structure paths. It can be 
easily implemented in existing codes, without searching to get cross-terms, 
reducing the numerical effort to adapt our method in these codes. 

Finally, as in the one-band case, it is necessary to find the best Slater 
determinant leading to optimized effective levels. One can easily show that 
expression (21) can be generalized for orbital degeneracy: 

2e . dln(Vq—) 

OTliaa 

Again, as in the one-band case, it is necessary to perform self-consistent 
calculations (for a given previously converged {po(T)} set) till the overall 
convergence is reached i.e. the on-site effective levels as well as the hoppings 
are converged. Once achieved, the effective Hamiltonian H e ff can be used 
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to a quasiparticles description of the system as proposed by Vollhardt [15]. 
In fact, he has shown that the Gutzwiller method is a natural frame to 
obtain the parameters of the phenomenological Landau theory of Fermi 
liquids. One can also expect finite temperature extension of our method in 
analogy with [47] and references therein. 



2.5. AB-INITIO APPROACH: ALTERNATIVE TO THE LDA+J7 METHOD 



We present now how to implement such an approach in an ab-initio calcula- 
tion of solids. The linearized muffin-tin orbital in the atomic sphere approx- 
imation (LMTO-ASA) is widely used and peculiarly suited for our purpose 
as the basis set has a local representation. Other ab-initio approaches could 
be used, and if they have not this local property, one could transform the 
basis into a Wannier representation. The LMTO method is well described 
elsewhere [29, 30] and we would like to remind here only the main results 
which are usefull for this paper. In the frame of DFT-LDA band structure 
calculations, the LMTO method is based on some approximations. The 
space is divided in atomic spheres where the potential is spherically sym- 
metric and interstitial region where it is flat ("Muffin Tin" potential). In 
the Atomic Sphere Approximation (A.S.A.), the spheres radii are chosen so 
that the total volume of the spheres equals that of the solid. One makes a 
further approximation by supposing that the kinetic energy in the intersti- 
tial region is zero (without this non-essential assumption, Laplace equation, 
as used below, should be replaced by Helmoltz equation). In this region, 
the Schrddinger equation reduces to Laplace equation having regular and 
irregular solutions: Yi,(f)r and Y£(f)r _1 respectively. Here L = (£,m) 
represents the angular momentum index and Yl{t) the spherical harmon- 
ics in direction f = (9,<p). For the sphere centered at site R and in the 
momentum index £ (I = 0, 1, 2, 3), one finds the solution tpRtu of the radial 
Schroedinger equation for a given energy E u , usually taken at the center of 
gravity of the occupied part of the £-band and the energy derivative of <pr£ v 
noted <j>R£u- It can be shown that the corresponding orbitals <pr£ U an d ^Riv 
are orthogonal to each other and nearly orthogonal to the core levels. It is 
thus possible to build a basis set of orbitals xrl centered at sphere of site 
R in the following way. Outside the sphere, in the interstitial region xrl is 
proportional to the irregular solution Yi,(r)r -1 of Laplace equation and 
it is augmented (i.e. substituted according to Slater terminology) in its own 
sphere by a linear combination of <piuu and <piuu having logarithmic deriva- 
tive — I — 1 at the radius sr of the sphere so that the orbital is continuous 
and derivable at the sphere boundary. In any other sphere R', the irregular 
solution of Laplace equation can be expanded in term of regular solutions 
in that sphere: 
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and the orbital xrl should be augmented in sphere R' with the same 
expansion of linear combination of <Pr'£' v and ^>B!i'v having the logarithmic 
derivative £' at the radius sr> of sphere R'. In (37), a is a scale factor and 
Sr'L' rl are ^ ne so-called "structure constants" which depend only on the 
crystallographic structure of the material. In this basis set of the orbitals 
Xrl both Hamiltonian and Overlap matrices can be expressed in terms of 
S R , L , RL , and the potential parameters <£riu{sr), tpRtv{ s R) an d their loga- 
rithmic derivatives Dri v and Dr£ u at sphere boundary. Since the structure 
constants S R , L , RL , decreasing as r~ e ~ e _1 with distance, are very long 
ranged for s and p orbitals, it can be more convenient to change the basis 
set so that the Hamiltonian can have the Tight-Binding (TB) form or any 
desired properties (like the orthogonality of overlap). It can be achieved by 
adding to the regular solution of Laplace equation an amount of the irreg- 
ular solution for a given angular momentum. It is possible to choose this 
amount Qi so that the transformed structure constants S can be screened 
with a short-range dependence with the distance or so that the orbitals of 
the transformed basis set are orthogonal (the so-called TB or most local- 
ized and orthogonal representations, respectively). With appropriate choice 
for Qg, the transformed structure constant matrix obeys to the following 
equation: 

S = S°(l - QeS )- 1 (38) 
Matrix elements fo the Hamiltonian can be written as: 

Hrl,r'L' = Crl<>rl,r'L' + Ar^Srl,r'u (39) 

which is limited to first order in {E — E v ) in the TB representation, whereas 
it is valid up to second order in the orthogonal representation (and it is even 
possible to add third order correction). C R l determines the middle of the 
band " RL" and A R l its width and the strength of hybridization. These 
parameters are expressed in terms of the 4 potential parameters: (Priu(sr), 
( PiUv(sr), Dri u and Dri v . It should be stressed that hybridization between 
bands of different angular moments is due to the matrix elements Srl^r'L' 
which couples i?L-states to R'L' ones. When these matrix elements are set 
equal to zero for £ / £', one obtains bands having pure £ character. This 
approximation was suggested in the standard (unscreened) representation 
and the resulting bands were called "canonical" bands. In that case it would 
be quite easy to apply single band Gutzwiller method (one equation per £ 
symmetry) without the need of the previous fully hybridized generalization. 
However, as we want to treat realistic bands, we do not use the canonical 
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bands in this paper. We use a scalar relativistic LMTO-ASA code neglect- 
ing spin-orbit coupling, with the so-called "combined corrections" which 
correct the ASA. The density functional formalism, before the Gutzwiller 
correction explained below, was treated within the LDA with the exchange 
and correlation potential of von Barth and Hedin [31]. 

The Hamiltonian of valence electrons (39), in the so-called orthogonal 
representation (or in the most localized representation, neglecting orbital 
overlap) can be mapped on a tight-binding form Hamiltonian 

(40) 

i^jafia tour 

The hoppings and on-site energies are directly outputs of the ab-initio cal- 
culation, as explained in Ref. [30] and by identification of expression (40) 
with (39), making the correspondance: R — > i, L — > a. This opens the 
possibility of treating our approach from first principle level, without any 
adjustable parameters, except the interactions U. They could be however 
also computed from constrained LDA calculations but in the following we 
rather treat them as free parameters. As in (32), (40) describes a full spd 
(and possibly / as in application of this method to Plutonium, see next 
section) basis. The terms ej QCT account for different on-site energies for or- 
bitals with different angular momentum or lying on inequivalent sites with 
possible crystal field splitting between orbitals of same angular momen- 
tum, but belonging to different irreducible group representations: it is due 
to the on-site contribution of Srl,rl that arises in the TB or in the nearly 
orthogonal representation. 

In the spirit of the Anderson model, we separate electrons into two 
subsystems: delocalized electrons for which the LDA is assumed to give 
reasonable results and localized electrons for which it is well known that 
the LDA can lead to unphysical results. To treat these states in a better 
way, and to avoid double counting, we exclude the interaction between 
localized electrons (/ or d) already taken into account in an average way 
in the LDA-on-site energy 

<L, = & ~ U{n f - \) (41) 

where U is a proper combination of direct and exchange Coulomb integral 
giving a true one-electron Hamiltonian Hq. rif is the average number of / 
(or d) electrons given by the LDA calculation. We then re-add an inter- 
action part .ffint " a la Hubbard" for the localized electrons and the full 
Hamiltonian H' = Hq + H mt is treated within the previously described 
multiband Gutzwiller approach. In fact this starting Hamiltonian H' is the 
same one used in the so-called "LDA+?7" method [8], the difference being in 
the way the interaction part is treated. In the LDA+J7 method, it is treated 
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in a mean field Hartree-Fock like approach, which can be questionable in 
case of strong correlations. It is however a suitable way of introducing an 
orbital-dependent potential which is absent in DFT formalism. In our ap- 
proach the correlation is treated exactly, within the approximation of the 
Gutzwiller ansatz. Note that our method also contains an orbital-depend 
potential through the renormalization of levels (36). A detailed study of 
the involved derivative indeed reveals an orbital filing dependence when 
rewriting this formula: 



Similarly to what happens in the LDA+C7 method, one sees the tendency 
of lowering for levels with occupation greater than one half, and a rising 
upwards for the less than half filled ones (the partial kinetic energy ei aa 
is always negative, it would be zero for a filled band). The difference with 
LDA+f/" method is the partial kinetic energy prefactor (instead of U) and 
other terms that come from the derivative of the empty and single occupied 
configurations. 

The starting Hamiltonian H 1 has been also used to make a link be- 
tween ab-initio LMTO band structure calculation and a DMFT treatment 
of correlations for the studies of LaTiC>3 [9] and Plutonium [10]. This last 
approach, assuming infinite dimension, goes beyond our approach. We only 
expect to be able to describe the coherent part of the spectrum, whereas 
the incoherent part leading to lower and upper Hubbard subbands are not 
accessible in our model, however as already stressed, variationally based. 

The practical scheme we proposed to perform our ab-initio Gutzwiller 
approach is the following one. First, we perform a LDA ab-initio LMTO- 
ASA calculation of the solid in a given crystal structure. This calculation 
provides the core and the valence (band) electrons contribution total energy, 
as well as occupations and partial kinetic energies for valence orbitals. From 
these ingredients, and for a given model interaction Hamiltonian, it is possi- 
ble to evaluate the variational Gutzwiller energy, which will be minimized, 
providing an optimized set of variational configurational probabilities. Then 
the on-site levels are varied according to the prescription renormalization 
of levels (36) as well as the adjunction of q- factors (27) which modified 
hoppings. New partial kinetic energies and occupations are recalculated 
from the modified Hamiltonian H', until self-consistency is achieved. At 
the end of procedure, the total energy, sum of the core and band energies, 
is calculated, leading to the properties of the ground state. One can then 
change, for example, the volume and redo the whole loop to obtained the 
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equilibrium properties and the most favorable atomic configurations in the 
solid. Close to the Fermi level we can also obtain an approximate ab-initio 
description of quasiparticles spectrum which enables comparison with spec- 
troscopy experiments for moderately correlated electron systems. That way, 
one has an ab-initio method which is multi-configurational and variational. 

In this first application, we made a slight simplification, with respect 
to the general process described in last paragraph, for the band calcula- 
tion part: starting from converged LMTO potential parameters, we build 
up a first order Hamiltonian in TB representation with neglect of overlap 
matrix (i.e. equal to unity) as explained in [32] and references therein. As 
our scheme only reorganizes valence (band) electrons, we make a frequently 
used frozen core approximation assuming that the core energy remains un- 
affected by this reorganization and we will now concentrate on the band 
energies. It is well-known that this first order Hamiltonian is accurate close 
to E v , i.e. close to the center of gravity of the occupied part of the bands. 
Far from it, it has the effect of a slight reduction of the bandwidth, but 
we have verified that it has a negligible effect on integrated quantities: for 
example, before the Gutzwiller process is switched on, we have checked that 
the band energies calculated from the third order Hamiltonian and from the 
first order one with the recursion process described below, are in excellent 
agreement. We used a full spdf basis set with hoppings up to second near- 
est neighbors. For all 7 inequivalent orbitals, in cubic environment from the 
overall 16 orbitals, we performed a real space recursion procedure [33] to 
get the partial projected densities of states (DOS) from which all needed 
quantities, like occupancies or band energies, can be calculated. These par- 
tial DOS are obtained from the imaginary part of diagonal elements of a 
Green function, which are developed in a continued fraction expansion up 
to a given level. This level is chosen so that a convergence criterium is 
reached, i.e. adding one more level does not affect the result. Practically we 
took 40 steps of recursion. Various terminators (the well-known square root 
terminator, or more elaborated ones in presence of gaps [34]) are then used 
to close the continued fraction expansion. A full self-consistent approach 
within the Gutzwiller loop, using third order Hamiltonians and including 
spin-orbit coupling, is still in preparation, and some intermediary results 
will be given in next section. 

3. Application to Plutonium 

We now give a simple application of the present method to Plutonium 
which is a good test case. Pu lies between light actinides with itinerant 
5/ electrons and heavy actinides with localized 5/ electrons. The compe- 
tition between these two electronic regimes in Pu is responsible for a lot 
of unusual properties as large values of the linear term in the specific heat 
coefficient and of the electrical resistivity or a very complex phase diagram. 
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The ground state a phase (monoclinic with 16 atoms by cell) is known 
to be well described by ab-initio DFT-LDA calculations, whereas for the 
high temperature 5 phase (fee) , the calculated equilibrium volume is of the 
order of 30 percent smaller than the experimental one. It is very important 
to reproduce the properties of Pu to take into account the delicate bal- 
ance between the itinerancy of the / electrons and the large intra atomic 
Coulomb interaction. This requires a much more complicated theory for the 
electrons than the LDA which is like a mean-field treatment of the correla- 
tions [35]. Recently several attempts to go beyond LDA have given a new 
understanding of the a- 5 transition. In the LDA+C/ method [36] an orbital- 
dependent correction, treated in the mean-field approximation is added to 
the LDA functional. These calculations have showed how the equilibrium 
volume is improved in comparison to previous results using LDA, and how 
an augmentation of the orbital moment is observed following Hund's rules, 
reducing the total magnetic moment in agreement with experiments. Going 
a step beyond the LDA+J7, Savrasov et al [37] have used an implementa- 
tion of DMFT. The LDA+[7 can be viewed as the static approximation of 
the DMFT. With this dynamical treatment of the /-electrons they have 
recovered the experimental equilibrium volume of e)-Pu, the photoemission 
peak at the Fermi level and given an understanding picture of the transi- 
tion between a and 5 phases. Different approaches, using the spin-polarized 
generalized gradient approximation (GGA) and antiferromagnetic configu- 
rations [38-40] have well reproduced the ground state properties of 5-Pu. 
All these works show how a spin/orbital polarization is crucial to describe 
the 5-phase. In the Gutzwiller method, the correlations, via the (/-factors, 
are supposed to reduce the hoppings, and so to weaken the covalency char- 
acter of the bonding, and consequently the attraction between atoms. Thus 
we expect to increase the interatomic distance, leading to a greater equi- 
librium volume. Of course, the same approach has to be performed for a 
and 5 phase. 

An extra difficulty arises from the Atomic Sphere Approximation (ASA) 
of the LMTO method: the atomic potential, inside an atomic "muffin-tin" 
sphere, is spherized, or equivalently, the true "full" potential is approxi- 
mated by its first £ = component. This approximation greatly simplifies 
the calculation, as the wave function basis in a sphere, used to build the 
LMTO set, can be factorized in a product of a radial wave function and a 
spherical harmonics as explained above. It presents however the shortcom- 
ings that, it is not a "full potential" approach and forbids to change the 
symmetry when making comparison between structures. We overcome this 
difficulty here by performing the calculation in a fee structure browsing dif- 
ferent volume: it is correct for the 5 phase, but the a phase will be replaced 
by a "pseudo"-a phase, in a fee structure, having however the same density 
than the experimental one. 

The valence states taken into account in the LMTO part were the 7s, Qp, 
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6d, 5/ of Pu with 16 fully hybridized orbitals per site, the remaining orbitals 
being treated as core states. In this first approach, as we concentrate more 
on the correlation effects, we neglect, however important for this heavy 
element, the spin-orbit coupling. The crystal field splitting on / orbitals 
(and other ones), is directly accounted by the LMTO method, lifting the 
/ degeneracy in the 6- fold (including spin) T\, the 6- fold T2 and 2-fold A2 
symmetries. Finally, the interaction H int , added to H , is simply given by 
the same local term between electrons on different /-orbitals 



neglecting any exchange term as done in [36,37]. In this simplified paramag- 
netic version, the number of inequivalent atomic configurations, necessary 
to perform the Gutzwiller part, reduces to 14 because all atomic configu- 
rations having the same electronic occupancy are equivalent in this model. 
Similarly, we took an average occupation per / orbital in the expression of 
g-factors, leading to a single q for all / orbitals, regardless to crystal field 
splitting. It was, however, included for the on-site levels renormalization, 
since the partial kinetic energy and occupations are not exactly equal for 
different symmetries. We have nevertheless checked this assumption by per- 
forming a much heavy calculation, including 3 different g's, one per crystal 
symmetry with 7x7x3 = 147 variational parameters: the final result was 
not sensitive to this detail. It reflects the small /-crystal field splitting in 
Plutonium, producing very similar occupations and partial kinetic energies. 

The Coulomb interaction U could be also provided by constrained LDA 
calculations. In that sense, it would not be an adjustable parameter. How- 
ever, we did not recalculate its value and took it from literature, close to 
0.3Ry, as in the LDA+DMFT calculation of Savrasov et al. [10], or as in the 
LDA+[7 calculation of Bouchet et al. [36]. An improved version of calcula- 
tion, including exchange interaction, as in the degenerate Hubbard model, 
with one (/-factor per symmetry, will be used in a forthcoming paper, in 
which we will investigate also ferromagnetic and antiferromagnetic ground 
states. In this work we just want to appreciate the effect of our method and 
of the Gutzwiller approximation on a simple case, where there exists known 
results with other methods. 

The total energy versus volume for fcc-Pu and different values of the 
interaction U is presented in Fig. 4. The curve U = corresponds to a 
LDA calculation. As previously found in several works the minimum of this 
curve is very low (~ 7.70 ua ) compared to the experimental value of the 
5 phase (8.60 ua) and closer to the a phase value (8.0 ua). In fact there is 
no sign of the correlated 5 phase in the U = calculation. As we turn on 
the correlations, a new feature appears in the curves, almost instantly. We 
observe a new energy minimum close to the experimental volume of the 5 
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phase. Moreover the first minimum increases to approach the value of the 
experimental a volume, showing that correlations are already important to 
reproduce the properties of this phase. For a value of U close to 0.3 Ry, the 
two minimums correspond to the experimental values of a and S-Pu. This 
double- well feature of the total energy curve of Pu was previously discovered 
by Savrasov et al [10], using a DMFT approach. In our calculations the 
first minimum is the lower one, since the a phase is the ground state for 
Pu, and we haven't added any temperature effect in our calculations. As 
U increases we see a tendency of the two minimums to be closer. In fact 
the energies of the two phases are very similar and a small perturbation, 
for example the temperature, can be sufficient for the phase transition. Of 
course the model studied in this work is still very simple and we don't want 
to conclude too far but we think that it already contains the key ingredients 
(competition between localization and derealization, atom-like or bands- 
like descriptions) to reproduce the main characteristics of Plutonium phase 
diagram. Due to the roughness of our first approach, the (rather) good 
agreement for the equilibrium properties, may be incidental or due to some 
compensation effect, and the disagreement with other aspects (like bulk 
modulus, see below) is not surprising. Indeed, it is well known the the spin- 
orbit coupling is a key ingredient for this element: the splitting between 5/2 
and 7/2 states could give significant differences in occupation and kinetic 
energies. One may expect then a difference between g 5 / 2 and q-j/2, and 
obtain localized and less localized behaviors as suggested by Penicaud [41] 
who proposed to split / states between localized and more delocalized ones 
to explain the properties of Plutonium. The freezing of /-states to similar 
occupation in our present calculation could be responsible for the high value 
of the bulk modulus (637 GPa) we get, in contrast with the experimental 
value of 30 GPa [42]. Primary result with an improved version involving 
third order LMTO Hamiltonian full self-consistent computation, neglecting 
yet spin-orbit coupling, reduces this value to 196 GPa, which is slightly 
better than the LDA result of 214 GPa [43]. 

This ab-initio Gutzwiller approach is able to handle correctly the cor- 
relation aspects without loosing the ab-initio adjustable parameters free 
aspect of the more familiar DFT-LDA, and that way, corrects the defi- 
ciency of this method. It gives similar results to the methods that account 
for many-body effects like the LDA+DMFT of Ref. [10] from the ab-initio 
levels or that can have an orbital dependent potential like in the LDA+C/ 
calculation of Ref. [36], which is impossible to DFT-LDA approach. On 
another hand, we stress again that our approach is clearly variational, and 
is able to provide an approximate ground state in contrast with those of 
Refs. [10] and [36]. 

The effective optimized Hamiltonian H', was used to compute quasipar- 
ticles density of states, in the vicinity of Fermi energy. The result, shown 
on Fig. 5, is restricted to an energy window of 2eV on both sides of Fermi 
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Figure 4- Total energy of fcc-Pu versus volume for different values of the interaction U. 
U = corresponds to a LDA calculation. 
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Figure 5. Quasiparticles density of states obtained from Gutzwiller method for Pluto- 
nium in S phase. 



level: further, the spectrum would be in the region of the Hubbard sub- 
bands, which are beyond the scope of our approach (it should be necessary 
to include fluctuations to get the incoherent part of the spectrum). It is 
to be stressed that our result compares well, in the presented region, with 
the more elaborated LDA+DMFT result of Ref. [10]. Both results are also 
in good agreement with the photoemission experiments of Arko et al. [44]. 
The peak at the Fermi level, which results mainly from the reduction of 
hoppings due to (/-factors, associated with a small shift of Fermi level due 
to the renormalization of levels (36), has the consequence of a significant 
improvement of the electronic specific heat contribution, multiplied by a 
factor of 2 with respect to the LDA value. Our result, of the order of 13 
mJ K -2 mollis however yet far from the experimental value of 64 mJ K~ 2 
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mol -1 found by Lashley et al [45]. The (/-factor, responsible for this increase 
of the DOS at the Fermi level, is of the order of 0.8 in the equilibrium 8 
phase. This moderate renormalization is due to hybridization of the / states 
with the low lying p-states which lead to a significant partial / kinetic en- 
ergy, greater than what it would be considering only the narrow group of 
predominant /-character states. With a naive single /-band argument, if 
we had used for example canonical (i.e. unhybridized) bands, one obtains a 
much reduced g-factor close to .3 [46]. It can be shown that the g-factor is 
the spectral weight of the quasiparticles. A moderate value, as obtained by 
our realistic calculation (i.e. fully hybridized bands), means a rather high 
weight: for independent particles it would be equal to one. It validates a 
quasiparticles picture description, allowing a posteriori comparison we did 
with spectroscopy experiments. At the volume of the "pseudo" a phase, 
this q factor reduces to .9, indicating that the electrons are less correlated 
in this phase, which can explain the relative success of its description by 
LDA calculation. 

4. Conclusion 

To conclude, we have generalized the density matrix approach to Gutzwiller 
method for the degenerate Hubbard Hamiltonian. We have shown that 
we can express the total energy in the Gutzwiller state in terms of the 
different probabilities of configurations. Moreover to apply the method to 
cases of physical interest we have developed this method for inequivalent 
sites and for different orbital symmetries. In this way we have given the 
expression of the different q factors which renormalize the hopping terms 
and an expression to renormalize the on-site energies in the Gutzwiller 
context. This method is limited to ground state properties but can extend 
to finite temperature and low-frequency excitations in analogy with the 
work of Gebhard [47]. Of course, as a quasiparticle approach, this method 
is limited to cases where the Fermi-liquid theory is valid, i.e. close to Fermi 
energy. Thereafter we have have described a simple implementation of our 
method in a ab-initio calculation as the LMTO method. To give an example, 
we have applied this technique to the particular case of Pu in fee structure. 
In despite of the simplicity of our model, we were able to extract interesting 
results such as the double- well feature in the energy- volume curve and more 
generally improve the LDA results. Our results compare well with previous 
works. 
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